Feature Selection in Machine Learning (Breast Cancer Datasets)

번역에 앞서.

이번 포스트는 최근 관심을 가지고 있는 Feature selection과 관련된 내용을 담고 있다. 완벽한 번역보다는 아티클을 보고 정리한 내용들을 거칠게 정리하려고 한다. 자세한 내용은 원문 포스트 참고하길 바란다.

Introduction

머신러닝은 예측 모델을 생성하기 위해서 변수(features, 또는 variables, attributes)를 사용한다. 높은 정확도를 얻기 위해서는 적절한 변수들의 조합을 사용하는 것이 필수적이다. (불명확한) 변수들을 사용하면 모델이 과적합되는 문제가 발생하기 때문에, 일반적으로 예측하고자 하는 반응변수와 가장 연관성이 높은 변수들을 찾아서 모델에 사용하길 원한다. 가능한한 적은 수의 변수를 사용하면 모델의 복잡도(complexity)를 낮출 수 있다. 다시 말해서 모델을 실행하는 시간과 컴퓨팅 파워 요구량을 줄이고, 이해하기도 쉬워진다는 의미다.

모델에 대해서 각각의 변수가 기여하는 정도를 측정하고, 변수의 수를 정하는 방법은 여러 가지가 있다. 본 아티클에서는 랜덤 포레스트 모델에 대해서 다음의 기법을 사용한다.

  • 상관계수 (Correlation),
  • 재귀 변수 제거법 (Recursive Feature Elimination),
  • 유전 알고리즘 (Genetic Algorithm, GA)

추가적으로, 데이터의 특성이 다를 때 위의 변수 선택 기법이 어떤 영향을 미치는지 확인하기 위해서 세 가지의 유방암 관련 데이터셋을 사용하였다. 하나는 매우 적은 수의 변수를 가지고 있으며, 다른 두 개의 데이터는 매우 큰 데이터지만, PCA를 이용해서 결과 클러스터를 다르게 하였다.

상관계수 기법, RFE, GA를 비교한 결과에 따르면, 랜덤 포레스트에 대해서 다음의 결과를 얻었다.

  • 상관계수가 높은 변수를 제거하는 것이 일반적으로 적절한 기법은 아니다.
  • GA가 본 예제에서는 최상의 결과를 보여줬지만, 다양한 변수가 있는 사례에서는 실용적이지 못하다. 적절한 세대(generation)까지 모델을 실행하는 데 오랜 시간이 걸리기 때문이다.
  • 시작부터 좋은 분류 결과가 나오기 힘든 데이터는 변수 선택에서 큰 효과를 보기 어렵다.

이 결론들은 물론 모든 데이터에 대해 일반화할 수는 없다. 위 기법들 이외에도 다양한 변수 선택 기법이 있으며, 굉장히 제한된 데이터셋에 대해서만 비교분석 하였으며, 랜덤 포레스트 모델에 대한 영향력만 분석했다. 하지만 작은 예제로도 서로 다른 변수와 파라미터가 예측 결과에 어떻게 영향을 미치는지 충분히 보여주고 있다는 점이 중요하다. 머신러닝에서는 이른바 “만능(one size fits all)”은 존재하지 않는다. 데이터를 유심히 살펴보는 것은 항상 가치있는 일이고, 모델이나 알고리즘에 대해서 생각하기 전에 데이터의 특징들에 익숙해지는 것이 중요하다. 주어진 데이터에 대해서 무언가 감을 잡았다면, 서로 다른 변수 선택 기법(또는 생성한 변수들), 모델, 파라미터들을 비교해보는 데에 시간을 투자하고, 마지막으로 다양한 머신러닝 알고리즘을 비교해보는 것이 큰 차이를 만들어낼 수 있다.

Breast Cancer Wisconsin (Diagnostic) Dataset

변수 선택 기법을 비교하기 위해 사용할 데이터는 Breast Cancer Wisconsin (Diagnostic) 데이터셋이다.

W.N. Street, W.H. Wolberg and O.L. Mangasarian. Nuclear feature extraction for breast tumor diagnosis. IS&T/SPIE 1993 International Symposium on Electronic Imaging: Science and Technology, volume 1905, pages 861-870, San Jose, CA, 1993.

O.L. Mangasarian, W.N. Street and W.H. Wolberg. Breast cancer diagnosis and prognosis via linear programming. Operations Research, 43(4), pages 570-577, July-August 1995.

W.H. Wolberg, W.N. Street, and O.L. Mangasarian. Machine learning techniques to diagnose breast cancer from fine-needle aspirates. Cancer Letters 77 (1994) 163-171.

W.H. Wolberg, W.N. Street, and O.L. Mangasarian. Image analysis and machine learning applied to breast cancer diagnosis and prognosis. Analytical and Quantitative Cytology and Histology, Vol. 17 No. 2, pages 77-87, April 1995.

W.H. Wolberg, W.N. Street, D.M. Heisey, and O.L. Mangasarian. Computerized breast cancer diagnosis and prognosis from fine needle aspirates. Archives of Surgery 1995;130:511-516.

W.H. Wolberg, W.N. Street, D.M. Heisey, and O.L. Mangasarian. Computer-derived nuclear features distinguish malignant from benign breast cytology. Human Pathology, 26:792–796, 1995.

데이터는 UC Irvine Machine Learning Repository에서 다운로드했다. 데이터셋의 변수들은 세포핵의 특성을 담고 있으며, 유방 세포 덩어리 세침흡인검사의 이미지 분석으로 생성되었다.

총 세 개의 데이터셋을 포함하고 있다. 첫번째 데이터셋은 아홉 개의 변수만 사용하고 있는 작은 데이터이며, 다른 두 개의 데이터는 각각 30개와 33개의 변수를 포함하고 있다. 두 데이터는 PCA로 생성되는 클러스터가 서로 다르다. 서로 다른 특성을 가지고 있는 데이터셋을 이용해서 서로 다른 변수 선택 기법의 효과를 살펴보고자 한다.

Breast Cancer Dataset 1

반응변수의 클래스는 다음과 같다.

  • Malignant (악성)
  • Benign (양성).

세포 특징에 관한 표현형은 다음과 같다.

  • Sample ID (code number)
  • Clump thickness
  • Uniformity of cell size
  • Uniformity of cell shape
  • Marginal adhesion
  • Single epithelial cell size
  • Number of bare nuclei
  • Bland chromatin
  • Number of normal nuclei
  • Mitosis
  • Classes, i.e. diagnosis

결측값들은 mice 패키지를 이용해서 처리한다.

bc_data <- read.table("data/breast/breast-cancer-wisconsin.data.txt",
                      header = FALSE, sep = ",")
colnames(bc_data) <- c("sample_code_number",
                       "clump_thickness",
                       "uniformity_of_cell_size",
                       "uniformity_of_cell_shape",
                       "marginal_adhesion", 
                       "single_epithelial_cell_size", 
                       "bare_nuclei", 
                       "bland_chromatin", 
                       "normal_nucleoli",
                       "mitosis",
                       "classes"
)
bc_data$classes <- ifelse(bc_data$classes == "2", "benign",
                          ifelse(bc_data$classes == "4", "malignant", NA))
bc_data[bc_data == "?"] <- NA

# how many NAs are in the data
length(which(is.na(bc_data)))
[1] 16
# impute missing data
library(mice)

bc_data[,2:10] <- apply(bc_data[, 2:10], 2, function(x) as.numeric(as.character(x)))
dataset_impute <- mice(bc_data[, 2:10],  print = FALSE)
bc_data <- cbind(bc_data[, 11, drop = FALSE], mice::complete(dataset_impute, 1))

bc_data$classes <- as.factor(bc_data$classes)

# how many benign and malignant cases are there?
summary(bc_data$classes)
   benign malignant 
      458       241 
str(bc_data)
'data.frame':   699 obs. of  10 variables:
 $ classes                    : Factor w/ 2 levels "benign","malignant": 1 1 1 1 1 2 1 1 1 1 ...
 $ clump_thickness            : num  5 5 3 6 4 8 1 2 2 4 ...
 $ uniformity_of_cell_size    : num  1 4 1 8 1 10 1 1 1 2 ...
 $ uniformity_of_cell_shape   : num  1 4 1 8 1 10 1 2 1 1 ...
 $ marginal_adhesion          : num  1 5 1 1 3 8 1 1 1 1 ...
 $ single_epithelial_cell_size: num  2 7 2 3 2 7 2 2 2 2 ...
 $ bare_nuclei                : num  1 10 2 4 1 10 10 1 1 1 ...
 $ bland_chromatin            : num  3 3 3 3 3 9 3 3 1 2 ...
 $ normal_nucleoli            : num  1 2 1 7 1 7 1 1 1 1 ...
 $ mitosis                    : num  1 1 1 1 1 1 1 1 5 1 ...
Breast Cancer Dataset 2

두번째 데이터셋의 반응변수의 클래스는 동일하다.

  • Malignant (악성)
  • Benign (양성).

데이터셋의 첫 두 칼럼은 다음과 같다.

  • Sample ID
  • Classes, i.e. diagnosis

각 세포핵에 대해서, 다음 10개의 특징이 측정되어 있다.

  • Radius (mean of all distances from the center to points on the perimeter)
  • Texture (standard deviation of gray-scale values)
  • Perimeter
  • Area
  • Smoothness (local variation in radius lengths)
  • Compactness (perimeter^2 / area - 1.0)
  • Concavity (severity of concave portions of the contour)
  • Concave points (number of concave portions of the contour)
  • Symmetry
  • Fractal dimension (“coastline approximation” - 1)

각각의 특징은 세 개의 기준을 측정되었다.

  • 평균
  • 표준편차
  • 가장 심각한 경우
bc_data_2 <- read.table("data/breast/wdbc.data.txt",
                        header = FALSE, sep = ",")

phenotypes <- rep(c("radius",
                    "texture",
                    "perimeter",
                    "area",
                    "smoothness",
                    "compactness",
                    "concavity",
                    "concave_points",
                    "symmetry",
                    "fractal_dimension"), 3)
types <- rep(c("mean", "se", "largest_worst"), each = 10)

colnames(bc_data_2) <- c("ID",
                         "diagnosis",
                         paste(phenotypes, types, sep = "_")
)

# how many NAs are in the data
length(which(is.na(bc_data_2)))
[1] 0
# how many benign and malignant cases are there?
summary(bc_data_2$diagnosis)
  B   M 
357 212 
str(bc_data_2)
'data.frame':   569 obs. of  32 variables:
 $ ID                             : int  842302 842517 84300903 84348301 84358402 843786 844359 84458202 844981 84501001 ...
 $ diagnosis                      : Factor w/ 2 levels "B","M": 2 2 2 2 2 2 2 2 2 2 ...
 $ radius_mean                    : num  18 20.6 19.7 11.4 20.3 ...
 $ texture_mean                   : num  10.4 17.8 21.2 20.4 14.3 ...
 $ perimeter_mean                 : num  122.8 132.9 130 77.6 135.1 ...
 $ area_mean                      : num  1001 1326 1203 386 1297 ...
 $ smoothness_mean                : num  0.1184 0.0847 0.1096 0.1425 0.1003 ...
 $ compactness_mean               : num  0.2776 0.0786 0.1599 0.2839 0.1328 ...
 $ concavity_mean                 : num  0.3001 0.0869 0.1974 0.2414 0.198 ...
 $ concave_points_mean            : num  0.1471 0.0702 0.1279 0.1052 0.1043 ...
 $ symmetry_mean                  : num  0.242 0.181 0.207 0.26 0.181 ...
 $ fractal_dimension_mean         : num  0.0787 0.0567 0.06 0.0974 0.0588 ...
 $ radius_se                      : num  1.095 0.543 0.746 0.496 0.757 ...
 $ texture_se                     : num  0.905 0.734 0.787 1.156 0.781 ...
 $ perimeter_se                   : num  8.59 3.4 4.58 3.44 5.44 ...
 $ area_se                        : num  153.4 74.1 94 27.2 94.4 ...
 $ smoothness_se                  : num  0.0064 0.00522 0.00615 0.00911 0.01149 ...
 $ compactness_se                 : num  0.049 0.0131 0.0401 0.0746 0.0246 ...
 $ concavity_se                   : num  0.0537 0.0186 0.0383 0.0566 0.0569 ...
 $ concave_points_se              : num  0.0159 0.0134 0.0206 0.0187 0.0188 ...
 $ symmetry_se                    : num  0.03 0.0139 0.0225 0.0596 0.0176 ...
 $ fractal_dimension_se           : num  0.00619 0.00353 0.00457 0.00921 0.00511 ...
 $ radius_largest_worst           : num  25.4 25 23.6 14.9 22.5 ...
 $ texture_largest_worst          : num  17.3 23.4 25.5 26.5 16.7 ...
 $ perimeter_largest_worst        : num  184.6 158.8 152.5 98.9 152.2 ...
 $ area_largest_worst             : num  2019 1956 1709 568 1575 ...
 $ smoothness_largest_worst       : num  0.162 0.124 0.144 0.21 0.137 ...
 $ compactness_largest_worst      : num  0.666 0.187 0.424 0.866 0.205 ...
 $ concavity_largest_worst        : num  0.712 0.242 0.45 0.687 0.4 ...
 $ concave_points_largest_worst   : num  0.265 0.186 0.243 0.258 0.163 ...
 $ symmetry_largest_worst         : num  0.46 0.275 0.361 0.664 0.236 ...
 $ fractal_dimension_largest_worst: num  0.1189 0.089 0.0876 0.173 0.0768 ...
Breast Cancer Dataset 3

세 번째 데이터셋의 반응변수의 클래스는 다음과 같다.

  • R: 재발
  • N: 재발하지 않음

데이터셋의 첫 두 칼럼은 다음과 같다.

  • Sample ID
  • Classes, i.e. outcome

각 세포핵에 대해서 두 번째 데이터셋과 동일한 특성과 측정기준으로 구성된 칼럼이 있으며, 추가적으로 다음의 칼럼이 있다.

  • Time (recurrence time if field 2 = R, disease-free time if field 2 = N)
  • Tumor size - diameter of the excised tumor in centimeters
  • Lymph node status - number of positive axillary lymph nodes observed at time of surgery
bc_data_3 <- read.table("data/breast/wpbc.data.txt",
                        header = FALSE, sep = ",")
colnames(bc_data_3) <- c("ID",
                         "outcome", 
                         "time",
                         paste(phenotypes, types, sep = "_"),
                         "tumor_size",
                         "lymph_node_status"
)

bc_data_3[bc_data_3 == "?"] <- NA

# how many NAs are in the data
length(which(is.na(bc_data_3)))
[1] 4
# impute missing data
library(mice)

bc_data_3[,3:35] <- apply(bc_data_3[,3:35], 2, function(x) as.numeric(as.character(x)))
dataset_impute <- mice(bc_data_3[,3:35],  print = FALSE)
bc_data_3 <- cbind(bc_data_3[, 2, drop = FALSE], mice::complete(dataset_impute, 1))

# how many recurring and non-recurring cases are there?
summary(bc_data_3$outcome)
  N   R 
151  47 
str(bc_data_3)
'data.frame':   198 obs. of  34 variables:
 $ outcome                        : Factor w/ 2 levels "N","R": 1 1 1 1 2 2 1 2 1 1 ...
 $ time                           : num  31 61 116 123 27 77 60 77 119 76 ...
 $ radius_mean                    : num  18 18 21.4 11.4 20.3 ...
 $ texture_mean                   : num  27.6 10.4 17.4 20.4 14.3 ...
 $ perimeter_mean                 : num  117.5 122.8 137.5 77.6 135.1 ...
 $ area_mean                      : num  1013 1001 1373 386 1297 ...
 $ smoothness_mean                : num  0.0949 0.1184 0.0884 0.1425 0.1003 ...
 $ compactness_mean               : num  0.104 0.278 0.119 0.284 0.133 ...
 $ concavity_mean                 : num  0.109 0.3 0.126 0.241 0.198 ...
 $ concave_points_mean            : num  0.0706 0.1471 0.0818 0.1052 0.1043 ...
 $ symmetry_mean                  : num  0.186 0.242 0.233 0.26 0.181 ...
 $ fractal_dimension_mean         : num  0.0633 0.0787 0.0601 0.0974 0.0588 ...
 $ radius_se                      : num  0.625 1.095 0.585 0.496 0.757 ...
 $ texture_se                     : num  1.89 0.905 0.611 1.156 0.781 ...
 $ perimeter_se                   : num  3.97 8.59 3.93 3.44 5.44 ...
 $ area_se                        : num  71.5 153.4 82.2 27.2 94.4 ...
 $ smoothness_se                  : num  0.00443 0.0064 0.00617 0.00911 0.01149 ...
 $ compactness_se                 : num  0.0142 0.049 0.0345 0.0746 0.0246 ...
 $ concavity_se                   : num  0.0323 0.0537 0.033 0.0566 0.0569 ...
 $ concave_points_se              : num  0.00985 0.01587 0.01805 0.01867 0.01885 ...
 $ symmetry_se                    : num  0.0169 0.03 0.0309 0.0596 0.0176 ...
 $ fractal_dimension_se           : num  0.00349 0.00619 0.00504 0.00921 0.00511 ...
 $ radius_largest_worst           : num  21.6 25.4 24.9 14.9 22.5 ...
 $ texture_largest_worst          : num  37.1 17.3 21 26.5 16.7 ...
 $ perimeter_largest_worst        : num  139.7 184.6 159.1 98.9 152.2 ...
 $ area_largest_worst             : num  1436 2019 1949 568 1575 ...
 $ smoothness_largest_worst       : num  0.119 0.162 0.119 0.21 0.137 ...
 $ compactness_largest_worst      : num  0.193 0.666 0.345 0.866 0.205 ...
 $ concavity_largest_worst        : num  0.314 0.712 0.341 0.687 0.4 ...
 $ concave_points_largest_worst   : num  0.117 0.265 0.203 0.258 0.163 ...
 $ symmetry_largest_worst         : num  0.268 0.46 0.433 0.664 0.236 ...
 $ fractal_dimension_largest_worst: num  0.0811 0.1189 0.0907 0.173 0.0768 ...
 $ tumor_size                     : num  5 3 2.5 2 3.5 2.5 1.5 4 2 6 ...
 $ lymph_node_status              : num  5 2 0 0 0 0 0 10 1 20 ...
Principal Component Analysis (PCA)

데이터셋의 차원과 분산에 대한 아이디어를 얻기 위해서, 데이터의 샘플과 변수에 대해서 PCA 플롯을 그리도록 한다. 첫 두 주성분(principal components, PCs)는 데이터 분산의 다수를 설명할 수 있는 두 주성분을 의미한다.

# function for PCA plotting
library(pcaGoPromoter)
library(ellipse)

pca_func <- function(data, groups, title, print_ellipse = TRUE) {
    
    # perform pca and extract scores
    pcaOutput <- pca(data, printDropped = FALSE, scale = TRUE, center = TRUE)
    pcaOutput2 <- as.data.frame(pcaOutput$scores)
    
    # define groups for plotting
    pcaOutput2$groups <- groups
    
    # when plotting samples calculate ellipses for plotting (when plotting features, there are no replicates)
    if (print_ellipse) {
        
        centroids <- aggregate(cbind(PC1, PC2) ~ groups, pcaOutput2, mean)
        conf.rgn  <- do.call(rbind, lapply(unique(pcaOutput2$groups), function(t)
            data.frame(groups = as.character(t),
                       ellipse(cov(pcaOutput2[pcaOutput2$groups == t, 1:2]),
                               centre = as.matrix(centroids[centroids$groups == t, 2:3]),
                               level = 0.95),
                       stringsAsFactors = FALSE)))
        
        plot <- ggplot(data = pcaOutput2, aes(x = PC1, y = PC2, group = groups, color = groups)) + 
            geom_polygon(data = conf.rgn, aes(fill = groups), alpha = 0.2) +
            geom_point(size = 2, alpha = 0.6) + 
            scale_color_brewer(palette = "Set1") +
            labs(title = title,
                 color = "",
                 fill = "",
                 x = paste0("PC1: ", round(pcaOutput$pov[1], digits = 2), "% variance"),
                 y = paste0("PC2: ", round(pcaOutput$pov[2], digits = 2), "% variance"))
        
    } else {
        
        # if there are fewer than 10 groups (e.g. the predictor classes) I want to have colors from RColorBrewer
        if (length(unique(pcaOutput2$groups)) <= 10) {
            
            plot <- ggplot(data = pcaOutput2, aes(x = PC1, y = PC2, group = groups, color = groups)) + 
                geom_point(size = 2, alpha = 0.6) + 
                scale_color_brewer(palette = "Set1") +
                labs(title = title,
                     color = "",
                     fill = "",
                     x = paste0("PC1: ", round(pcaOutput$pov[1], digits = 2), "% variance"),
                     y = paste0("PC2: ", round(pcaOutput$pov[2], digits = 2), "% variance"))
            
        } else {
            
            # otherwise use the default rainbow colors
            plot <- ggplot(data = pcaOutput2, aes(x = PC1, y = PC2, group = groups, color = groups)) + 
                geom_point(size = 2, alpha = 0.6) + 
                labs(title = title,
                     color = "",
                     fill = "",
                     x = paste0("PC1: ", round(pcaOutput$pov[1], digits = 2), "% variance"),
                     y = paste0("PC2: ", round(pcaOutput$pov[2], digits = 2), "% variance"))
            
        }
    }
    
    return(plot)
    
}

library(gridExtra)
library(grid)
Dataset 1
p1 <- pca_func(data = t(bc_data[, 2:10]),
               groups = as.character(bc_data$classes),
               title = "Breast cancer dataset 1: Samples"
)
p2 <- pca_func(data = bc_data[, 2:10],
               groups = as.character(colnames(bc_data[, 2:10])), 
               title = "Breast cancer dataset 1: Features", 
               print_ellipse = FALSE)
grid.arrange(p1, p2, ncol = 2)

h_1 <- hclust(dist(t(bc_data[, 2:10]),
                   method = "euclidean"),
              method = "complete")
plot(h_1)

library(tidyr)
bc_data_gather <- bc_data %>%
    gather(measure, value, clump_thickness:mitosis)

ggplot(data = bc_data_gather,
       aes(x = value, fill = classes, color = classes)) +
    geom_density(alpha = 0.3, size = 1) +
    geom_rug() +
    scale_fill_brewer(palette = "Set1") +
    scale_color_brewer(palette = "Set1") +
    facet_wrap( ~ measure, scales = "free_y", ncol = 3)

Dataset 2
p1 <- pca_func(data = t(bc_data_2[, 3:32]),
               groups = as.character(bc_data_2$diagnosis),
               title = "Breast cancer dataset 2: Samples")
p2 <- pca_func(data = bc_data_2[, 3:32],
               groups = as.character(colnames(bc_data_2[, 3:32])),
               title = "Breast cancer dataset 2: Features",
               print_ellipse = FALSE)
grid.arrange(p1, p2, ncol = 2, widths = c(0.4, 0.6))

h_2 <- hclust(dist(t(bc_data_2[, 3:32]),
                   method = "euclidean"),
              method = "complete")
plot(h_2)

bc_data_2_gather <- bc_data_2[, -1] %>%
    gather(measure, value, radius_mean:fractal_dimension_largest_worst)

ggplot(data = bc_data_2_gather, aes(x = value, fill = diagnosis, color = diagnosis)) +
    geom_density(alpha = 0.3, size = 1) +
    geom_rug() +
    scale_fill_brewer(palette = "Set1") +
    scale_color_brewer(palette = "Set1") +
    facet_wrap( ~ measure, scales = "free_y", ncol = 3)

Dataset 3
p1 <- pca_func(data = t(bc_data_3[, 2:34]),
               groups = as.character(bc_data_3$outcome), 
               title = "Breast cancer dataset 3: Samples")
p2 <- pca_func(data = bc_data_3[, 2:34], 
               groups = as.character(colnames(bc_data_3[, 2:34])),
               title = "Breast cancer dataset 3: Features", 
               print_ellipse = FALSE)
grid.arrange(p1, p2, ncol = 2, widths = c(0.4, 0.6))

h_3 <- hclust(dist(t(bc_data_3[,2:34]),
                   method = "euclidean"),
              method = "complete")
plot(h_3)

데이터셋 1과 데이터셋 2는 양성과 음성을 잘 분류한다. 또한 해당 변수들에 기반을 둔 모델은 클래스 예측 성능이 좋을 것으로 보인다. 하지만 데이터셋 3은 서로 다른 그룹으로 군집화하지 못하는데, 이는 해당 변수들을 사용했을 때 예측 성능이 떨어질 것으로 예상된다.

데이터셋 2와 데이터셋 3의 변수들은 잘 구별되게 군집화되지 않는다. 많은 변수들이 유사한 패턴을 보이기 때문인 것으로 보인다. 따라서 세 개의 데이터셋에 대해서 적절한 변수 부분집합을 고르는 것은 서로 다른 효과를 보일 것으로 보인다.

bc_data_3_gather <- bc_data_3 %>%
    gather(measure, value, time:lymph_node_status)

ggplot(data = bc_data_3_gather, aes(x = value, fill = outcome, color = outcome)) +
    geom_density(alpha = 0.3, size = 1) +
    geom_rug() +
    scale_fill_brewer(palette = "Set1") +
    scale_color_brewer(palette = "Set1") +
    facet_wrap( ~ measure, scales = "free_y", ncol = 3)

Feature importance

변수 각각의 중요성에 대한 정보를 얻기 위해서 caret 패키지를 사용하여 랜덤 포레스트에 대해 10 x 10 CV를 수행하였다. 모델링을 위한 변수 선택을 위해 변수 중요성이 필요했다면, 완전한 데이터셋이 아닌 트레이닝 데이터에 CV를 수행하여야 한다. 하지만 데이터에 전체에 대한 정보를 얻고 싶었기 때문에 전체를 사용하였다.

library(caret)
library(doMC)
registerDoMC(cores = 4)

# prepare training scheme
control <- trainControl(method = "repeatedcv", number = 10, repeats = 10)

feature_imp <- function(model, title) {
    
    # estimate variable importance
    importance <- varImp(model, scale = TRUE)
    
    # prepare dataframes for plotting
    importance_df_1 <- importance$importance
    importance_df_1$group <- rownames(importance_df_1)
    
    importance_df_2 <- importance_df_1
    importance_df_2$Overall <- 0
    
    importance_df <- rbind(importance_df_1, importance_df_2)
    
    plot <- ggplot() +
        geom_point(data = importance_df_1,
                   aes(x = Overall, y = group, color = group),
                   size = 2) +
        geom_path(data = importance_df,
                  aes(x = Overall, y = group,
                      color = group, group = group), 
                  size = 1) +
        theme(legend.position = "none") +
        labs(
            x = "Importance",
            y = "",
            title = title,
            subtitle = "Scaled feature importance",
            caption = "\nDetermined with Random Forest and
      repeated cross validation (10 repeats, 10 times)"
        )
    
    return(plot)
    
}
# train the model
set.seed(27)
imp_1 <- train(classes ~ .,
               data = bc_data,
               method = "rf",
               preProcess = c("scale", "center"),
               trControl = control)
p1 <- feature_imp(imp_1, title = "Breast cancer dataset 1")
set.seed(27)
imp_2 <- train(diagnosis ~ .,
               data = bc_data_2[, -1],
               method = "rf",
               preProcess = c("scale", "center"),
               trControl = control)
p2 <- feature_imp(imp_2, title = "Breast cancer dataset 2")
set.seed(27)
imp_3 <- train(outcome ~ .,
               data = bc_data_3,
               method = "rf",
               preProcess = c("scale", "center"),
               trControl = control)
p3 <- feature_imp(imp_3, title = "Breast cancer dataset 3")
grid.arrange(p1, p2, p3, ncol = 3, widths = c(0.3, 0.35, 0.35))